// ShekelFunction.cpp: implementation of the ShekelFunction class.
//
//////////////////////////////////////////////////////////////////////
#include <iomanip>
#include <fstream>
#include <string>

#include "shekel_function.h"
#include "DataHolder.h"



//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

ShekelFunction::ShekelFunction(int i_NInput)
{
  NInput = i_NInput;

  // Create 3 arrays
  newA(b, 10);
  newM(c, NInput, 10);
  newA(dDerivative, 10);

  newA(danswer, 10);
  newM(dDeriv,10,NInput);

  newA(mDerivative, NInput);

  Result.Size(NInput);

  knownMin.push_back(4.00075);
  knownMin.push_back(4.00059);
  knownMin.push_back(3.99966);
  knownMin.push_back(3.99951);
}


ShekelFunction::~ShekelFunction()
{
  deleteA(b);
  deleteM(c);
  deleteA(dDerivative);
  deleteA(mDerivative);
  
  deleteA(danswer);
  deleteM(dDeriv);
}


// takes a point and if it is beyong boundary move the point
// into the boundary box
void ShekelFunction::checkBoundary(vector<double> &Input)
{
  for (unsigned int i = 0; i < Input.size(); i++)
  {
    if ( Input[i] >= 10.0 )
    {
      Input[i] = fmod(Input[i],10.0);
    }
    else if ( Input[i] <= 0.0 )
    {
      Input[i] = fmod(Input[i],10.0)+10.0;
    }
  }
}

// takes a point and if it is beyong boundary move the point
// into the boundary box
void ShekelFunction::checkBoundary(double *Input)
{
  for (int i = 1; i <= NInput; i++)
  {
    // reflective boundary
    /*while( (Input[i] >= 10.0) || (Input[i] <= 0.0) )
    {
      if ( Input[i] >= 10.0 )
      {
        Input[i] = 10.0 - fabs(Input[i]-10.0);
      }
      else if ( Input[i] <= 0.0 )
      {
        Input[i] = 0.0 + fabs(0.0-Input[i]);
      }
    }*/

    // periodic boundary
    if ( Input[i] >= 10.0 )
    {
      Input[i] = fmod(Input[i],10.0);
    }
    else if ( Input[i] <= 0.0 )
    {
      Input[i] = 10.0 + fmod(Input[i],10.0);
    }

    // hard boundaries
    /*if ( Input[i] >= 10.0 )
    {
      Input[i] = 10.0;
    }
    else if ( Input[i] <= 0.0 )
    {
      Input[i] = 0.0;
    }*/
  }
}

//////////////////////////////////////////////////////////////////////
// functions
//////////////////////////////////////////////////////////////////////

double ShekelFunction::Value(double *Input)
{
  // first check that point is not beyond any boundaries
  checkBoundary(Input);

  // Initialise vars answer and danswer (both local)
  double answer = 0.0;
  //double *danswer;

  //newA(danswer, 10);

  if(NInput == 2)
  {
    c[1][1] = 4.0; c[1][2] = 1.0; c[1][3] = 8.0; c[1][4] = 6.0; c[1][5] = 3.0; c[1][6] = 2.0; c[1][7] = 5.0;
    c[1][8] = 8.0; c[1][9] = 6.0; c[1][10] = 7.0;
    c[2][1] = 4.0; c[2][2] = 1.0; c[2][3] = 8.0; c[2][4] = 6.0; c[2][5] = 7.0; c[2][6] = 9.0; c[2][7] = 5.0;
    c[2][8] = 1.0; c[2][9] = 2.0; c[2][10] = 3.6;
  } else {
    c[1][1] = 4.0; c[1][2] = 1.0; c[1][3] = 8.0; c[1][4] = 6.0; c[1][5] = 3.0; c[1][6] = 2.0; c[1][7] = 5.0;
    c[1][8] = 8.0; c[1][9] = 6.0; c[1][10] = 7.0;
    c[2][1] = 4.0; c[2][2] = 1.0; c[2][3] = 8.0; c[2][4] = 6.0; c[2][5] = 7.0; c[2][6] = 9.0; c[2][7] = 5.0;
    c[2][8] = 1.0; c[2][9] = 2.0; c[2][10] = 3.6;
    c[3][1] = 4.0; c[3][2] = 1.0; c[3][3] = 8.0; c[3][4] = 6.0; c[3][5] = 3.0; c[3][6] = 2.0; c[3][7] = 3.0;
    c[3][8] = 8.0; c[3][9] = 6.0; c[3][10] = 7.0;
    c[4][1] = 4.0; c[4][2] = 1.0; c[4][3] = 8.0; c[4][4] = 6.0; c[4][5] = 7.0; c[4][6] = 9.0; c[4][7] = 3.0;
    c[4][8] = 1.0; c[4][9] = 2.0; c[4][10] = 3.6;
  }

  b[1] = 0.1; b[2] = 0.2; b[3] = 0.2; b[4] = 0.4; b[5] = 0.4; b[6] = 0.6; b[7] = 0.3; b[8] = 0.7; b[9] = 0.5;
  b[10] = 0.5;

  for (int i = 1; i <= 10; i++)
  {
    danswer[i] = 0.0;
  }

  // Makes answer equal to it self added to
  // the sum of x^2/4000 for each dimension
  for (int n = 1; n <= 10; n++)
  {
    for (int m = 1; m <=NInput; m++)
    {
      danswer[n] += ((Input[m]-c[m][n])*(Input[m]-c[m][n]));
    }
    
    danswer[n] += b[n];
    answer -= 1.0/danswer[n];
  }

  Result.Input(answer, Input);

  nCall++;

  //deleteA(danswer);

  return (answer); 



}

// Calculates/builds the function value
double ShekelFunction::Value(vector<double> &Input)
{
  // first check that point is not beyong any boundaries
  checkBoundary(Input);

  // Initialise vars answer and danswer (both local)
  double answer = 0.0;
  //double *danswer;

  //newA(danswer, 10);

  if(NInput == 2)
  {
    c[1][1] = 4.0; c[1][2] = 1.0; c[1][3] = 8.0; c[1][4] = 6.0; c[1][5] = 3.0; c[1][6] = 2.0; c[1][7] = 5.0;
    c[1][8] = 8.0; c[1][9] = 6.0; c[1][10] = 7.0;
    c[2][1] = 4.0; c[2][2] = 1.0; c[2][3] = 8.0; c[2][4] = 6.0; c[2][5] = 7.0; c[2][6] = 9.0; c[2][7] = 5.0;
    c[2][8] = 1.0; c[2][9] = 2.0; c[2][10] = 3.6;
  } else {
    c[1][1] = 4.0; c[1][2] = 1.0; c[1][3] = 8.0; c[1][4] = 6.0; c[1][5] = 3.0; c[1][6] = 2.0; c[1][7] = 5.0;
    c[1][8] = 8.0; c[1][9] = 6.0; c[1][10] = 7.0;
    c[2][1] = 4.0; c[2][2] = 1.0; c[2][3] = 8.0; c[2][4] = 6.0; c[2][5] = 7.0; c[2][6] = 9.0; c[2][7] = 5.0;
    c[2][8] = 1.0; c[2][9] = 2.0; c[2][10] = 3.6;
    c[3][1] = 4.0; c[3][2] = 1.0; c[3][3] = 8.0; c[3][4] = 6.0; c[3][5] = 3.0; c[3][6] = 2.0; c[3][7] = 3.0;
    c[3][8] = 8.0; c[3][9] = 6.0; c[3][10] = 7.0;
    c[4][1] = 4.0; c[4][2] = 1.0; c[4][3] = 8.0; c[4][4] = 6.0; c[4][5] = 7.0; c[4][6] = 9.0; c[4][7] = 3.0;
    c[4][8] = 1.0; c[4][9] = 2.0; c[4][10] = 3.6;
  }

  b[1] = 0.1; b[2] = 0.2; b[3] = 0.2; b[4] = 0.4; b[5] = 0.4; b[6] = 0.6; b[7] = 0.3; b[8] = 0.7; b[9] = 0.5;
  b[10] = 0.5;

  // Constructing the equation
  for (int i = 1; i <= 10; i++)
  {
    danswer[i] = 0.0;
  }

  for (int n = 1; n <= 10; n++)
  {
    for (int m = 0; m <NInput; m++)
    {
      danswer[n] += ((Input[m]-c[m+1][n])*(Input[m]-c[m+1][n]));
    }
    
    danswer[n] += b[n];
    answer -= 1/danswer[n];
  }

  Result.Input(answer, Input);

  nCall++;

  //deleteA(danswer);

  return (answer); 

  

}

// Calculates/builds the derivatives
void ShekelFunction::BuildDerivative(double *Input, double *Derivative)
{
  //double **dDeriv;
  //newM(dDeriv,10,NInput);

  for (int i = 1; i <= 10; i++)
  {
    dDerivative[i]=0.0;
  }

  for (int j = 1; j <= NInput; j++)
  {
    Derivative[j]=0.0;
  }
  
  for (int k = 1; k <= 10; k++)
  {
    for (int l = 1; l <=NInput; l++)
    {
      dDerivative[k] += ((Input[l]-c[l][k])*(Input[l]-c[l][k]));
    }
    
    dDerivative[k] += b[k];
  }

  for (int m=1; m<=NInput; m++)
  {
    for(int n=1; n<=10; n++)
    {
      dDeriv[n][m] = (Input[m]-c[m][n])*(1/(dDerivative[n]*dDerivative[n]));
    }
  }

  for (int p=1; p<=NInput; p++)
  {
    for (int r=1; r<=10; r++)
    {
      Derivative[p] += (2*dDeriv[r][p]);
      //Derivative[p] += dDeriv[r][p];
    }
  }

  for(int mw=0; mw<=NInput; mw++)
  {
    mDerivative[mw] = Derivative[mw];
  }

//  deleteM(dDeriv);
 
}

// used to test the function is working
void ShekelFunction::test_function()
{
  ofstream out("output/test_of_shekel.txt");
  
  out << "NUMBER_OF_DIMENSIONS: " << NInput << endl 
    << "      Subscript[x,k]      |     Value    |         Derivative[k]" << endl << endl;

  int NRun = 3;
  double *mInput;
  double *mDeriv;
  newA(mInput, NInput);
  newA(mDeriv, NInput);

  for(int run=1; run<=NRun; run++)
  {
    for(int n=1; n<=NInput; n++)
    {
      mInput[n]=(10.0*ran1());
    }

    for(int k=1; k<=NInput; k++)
    {
      out << setw(11) << mInput[k];
    }
    out << setw(16) << Value(mInput);

    BuildDerivative(mInput, mDeriv);

    for(int l=1; l<=NInput; l++)
    {
      out << setw(16) << mDerivative[l];
    }

    out << "\n\n";
    
  }

  deleteA(mInput);
  deleteA(mDeriv);
  out.close();
}

// used to print the solution (this is the old version, the new version is below)
void ShekelFunction::print_sol(double *solution, double funcVal, char filename[])
{
  ofstream out(filename);	

  out << "Number of dimensions: " << NInput << endl;
  out << "Final Shekel value = " << funcVal << ";" << endl;
  out << "Solution coordinates:" << endl;

  for (int i = 1; i <= NInput; i++)
  {
    out << solution[i] << "  ";
  }
 
  out << endl;

  out.close();
}

// used to print the solution to a file
void ShekelFunction::print_sol(char filename[])
{
  ofstream out(filename);	

  // getBestValue() returns the best value stored in the tracker
  out << "Number of dimensions: " << NInput << endl;
  out << "Final Shekel value = " << Result.getBestValue() << endl;
  out << "Solution coordinates:" << endl;
  
  // prints out the solution coordinates to a file
  Result.printSolution(out);

  out.close();
}

// distance from know global minimum which is located at the origin
double ShekelFunction::distanceFromGlobalMinimum(double *solution)
{
  double toReturn = 0.0;

  for (int i = 1; i <= NInput; i++)
  {
    toReturn += (solution[i]-knownMin[i-1])*(solution[i]-knownMin[i-1]);
  }

  double toReturnS = sqrt(toReturn);
  return toReturnS;
}